# clean environment and load libraries
rm(list=ls())
library(readxl)
library(expm) 

# load wide dataset
d <- read_excel("/replication/data/wide_data.xlsx")
 

# Keep full dataset
d = d[d$jurisdiction1787!="NA",]
d = d[d$jurisdiction1352!="NA",]
d <- d[is.na(d$jurisdiction1787)==FALSE, ]

# Total transitions
table(d$jurisdiction1352,d$jurisdiction1787)

# Recode factor variables
d$j.1352 = factor(d$jurisdiction1352,
                  levels = c("Royal Jurisdiction",
                             "Seigneurial Jurisdiction",
                             "Ecclesiastical Jurisdiction",
                             "Elective Jurisdiction",
                             "Military Order",
                             "Mixed Jurisdiction"))

levels(d$j.1352) = c("R","S","E","B","O","M")

d$j.1787 = factor(d$jurisdiction1787,
                  levels = c("Royal Jurisdiction",
                             "Seigneurial Jurisdiction",
                             "Ecclesiastical Jurisdiction",
                             "Elective Jurisdiction",
                             "Military Order",
                             "Mixed Jurisdiction"))

levels(d$j.1787) = c("R","S","E","B","O","M")



# MAKE BARPLOTS
dev.off()
bar1 = cbind(d$jurisdiction1352,1352)
bar2 = cbind(d$jurisdiction1787,1787)
bar3 = rbind(bar1,bar2)
colnames(bar3) = c("Jurisdiction","Year")
bar3=as.data.frame(bar3)
bar3$Jurisdiction = factor(bar3$Jurisdiction,
                           levels = c("Royal Jurisdiction",
                                      "Seigneurial Jurisdiction",
                                      "Ecclesiastical Jurisdiction",
                                      "Elective Jurisdiction",
                                      "Military Order",
                                      "Mixed Jurisdiction"))

levels(bar3$Jurisdiction) = c("Royal","Seigneurial",
                              "Ecclesiastical","Behetria",
                              "Military Order","Mixed")

barplot(table(bar3$Year,bar3$Jurisdiction), 
        beside = TRUE,
        ylim = c(0,800))

legend("topright",
       c("1352","1787"),
       fill = c("black","grey"),
       cex = 1.3)



# TRANSITION MATRIX

total.trans = as.matrix(table(d$j.1352,d$j.1787))
total.trans

rowtotal = apply(total.trans, 1, sum)
rowtotal

transition.matrix = total.trans/rowtotal
transition.matrix

mat = matrix(NA, nrow = 6, ncol = 6)
mat[1, ] = transition.matrix[1, ]
mat[2, ] = transition.matrix[2, ]
mat[3, ] = transition.matrix[3, ]
mat[4, ] = transition.matrix[4, ]
mat[5, ] = transition.matrix[5, ]
mat[6, ] = transition.matrix[6, ]


########################################################
###              BOOTSTRAP STANDARD ERRORS             #
########################################################

list.mats = vector("list",10000)

set.seed(123)
for(i in 1:10000){

d.boot = d[sample(1:1624, 1624, replace = TRUE), ]  
total.trans = as.matrix(table(d.boot$j.1352,d.boot$j.1787))
total.trans

rowtotal = apply(total.trans, 1, sum)
rowtotal

transition.matrix2 = total.trans/rowtotal
transition.matrix2

mat = matrix(NA, nrow = 6, ncol = 6)
mat[1, ] = transition.matrix2[1, ]
mat[2, ] = transition.matrix2[2, ]
mat[3, ] = transition.matrix2[3, ]
mat[4, ] = transition.matrix2[4, ]
mat[5, ] = transition.matrix2[5, ]
mat[6, ] = transition.matrix2[6, ]

list.mats[[i]]=mat

}

# Compute bootstrap distribution

RR = c()
RS = c()
RE = c()
RB = c()
RO = c()
RM = c()

SR = c()
SS = c()
SE = c()
SB = c()
SO = c()
SM = c()

ER = c()
ES = c()
EE = c()
EB = c()
EO = c()
EM = c()

BR = c()
BS = c()
BE = c()
BB = c()
BO = c()
BM = c()

OR = c()
OS = c()
OE = c()
OB = c()
OO = c()
OM = c()

MR = c()
MS = c()
ME = c()
MB = c()
MO = c()
MM = c()



for(i in 1:10000){

 
 RR[i] = list.mats[[i]][1,1]
 RS[i] = list.mats[[i]][1,2]
 RE[i] = list.mats[[i]][1,3]
 RB[i] = list.mats[[i]][1,4]
 RO[i] = list.mats[[i]][1,5]
 RM[i] = list.mats[[i]][1,6]

 SR[i] = list.mats[[i]][2,1]
 SS[i] = list.mats[[i]][2,2]
 SE[i] = list.mats[[i]][2,3]
 SB[i] = list.mats[[i]][2,4]
 SO[i] = list.mats[[i]][2,5]
 SM[i] = list.mats[[i]][2,6]
 
 ER[i] = list.mats[[i]][1,1]
 ES[i] = list.mats[[i]][3,2]
 EE[i] = list.mats[[i]][3,3]
 EB[i] = list.mats[[i]][3,4]
 EO[i] = list.mats[[i]][3,5]
 EM[i] = list.mats[[i]][3,6]
 
 BR[i] = list.mats[[i]][4,1]
 BS[i] = list.mats[[i]][4,2]
 BE[i] = list.mats[[i]][4,3]
 BB[i] = list.mats[[i]][4,4]
 BO[i] = list.mats[[i]][4,5]
 BM[i] = list.mats[[i]][4,6]
 
 OR[i] = list.mats[[i]][5,1]
 OS[i] = list.mats[[i]][5,2]
 OE[i] = list.mats[[i]][5,3]
 OB[i] = list.mats[[i]][5,4]
 OO[i] = list.mats[[i]][5,5]
 OM[i] = list.mats[[i]][5,6]
 
 MR[i] = list.mats[[i]][6,1]
 MS[i] = list.mats[[i]][6,2]
 ME[i] = list.mats[[i]][6,3]
 MB[i] = list.mats[[i]][6,4]
 MO[i] = list.mats[[i]][6,5]
 MM[i] = list.mats[[i]][6,6]
 
}


# Bootstrap standard errors

m = matrix(NA, ncol = 6, nrow = 12)
m[1,] = round(transition.matrix[1, ], d=2)
m[2,1] = round(sd(RR), digits=2)
m[2,2] = round(sd(RS), digits=2)
m[2,3] = round(sd(RE), digits=2)
m[2,4] = round(sd(RB), digits=2)
m[2,5] = round(sd(RO), digits=2)
m[2,6] = round(sd(RM), digits=2)

m[3,] = round(transition.matrix[2, ], d=2)
m[4,1] = round(sd(SR), digits=2)
m[4,2] = round(sd(SS), digits=2)
m[4,3] = round(sd(SE), digits=2)
m[4,4] = round(sd(SB), digits=2)
m[4,5] = round(sd(SO), digits=2)
m[4,6] = round(sd(SM), digits=2)

m[5,] = round(transition.matrix[3, ], d=2)
m[6,1] = round(sd(ER), digits=2)
m[6,2] = round(sd(ES), digits=2)
m[6,3] = round(sd(EE), digits=2)
m[6,4] = round(sd(EB), digits=2)
m[6,5] = round(sd(EO), digits=2)
m[6,6] = round(sd(EM), digits=2)

m[7, ] = round(transition.matrix[4, ], d=2)
m[8,1] = round(sd(BR), digits=2)
m[8,2] = round(sd(BS), digits=2)
m[8,3] = round(sd(BE), digits=2)
m[8,4] = round(sd(BB), digits=2)
m[8,5] = round(sd(BO), digits=2)
m[8,6] = round(sd(BM), digits=2)

m[9, ] = round(transition.matrix[5, ], d=2)
m[10,1] = round(sd(OR), digits=2)
m[10,2] = round(sd(OS), digits=2)
m[10,3] = round(sd(OE), digits=2)
m[10,4] = round(sd(OB), digits=2)
m[10,5] = round(sd(OO), digits=2)
m[10,6] = round(sd(OM), digits=2)

m[11, ] = round(transition.matrix[6, ], d=2)
m[12,1] = round(sd(MR), digits=2)
m[12,2] = round(sd(MS), digits=2)
m[12,3] = round(sd(ME), digits=2)
m[12,4] = round(sd(MB), digits=2)
m[12,5] = round(sd(MO), digits=2)
m[12,6] = round(sd(MM), digits=2)

colnames(m) = c("R","S","E","B","O","M")

library(xtable)
xtable(m)

